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\\w problem of wave pr()pagati(ui in a simple viscoelastic string subjected to constant velocity 
transverse impact is reexamined. An error in an earlier solution by Smith (J. Res. NHS 70B, 257 
(1966)) is corrected and an alternate numerical scheme based on the method of characteristics but 
using an implicit formulation of finite-difference equations is presented. The constitutive ecjuation 
used is that of a linear viscoelastic model consisting of a spring and Maxwell element in parallel. Re- 
sults of an illustrative calculation are discussed. 
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1. Introduction 

The theoretical problem of wave propagation in a uniform string subjected to oblique impact 
is well understood if one assumes a dynamic stress-strain relation T= Tie) for the material, where 
T denotes the tensile forfce acting along the string, e, the Lagrangian strain, and there is no explicit 
dependence of T on the rate of strain. 

Several types of mathematical compHcation can enter into this theory. The first type is due 
to the nimlinearity of the function T{e). The theory predicts that: 

(a) if the curve Tie) is concave towards the e axis, strains of increasing magnitude are propa- 
gated with decreasing velocity c given by 



m de 

where m denotes the mass of the string per unit of unstrained length, and 

(b) if the curve T(e) contains a portion that is concave towards the T axis the strain distri- 
bution is discontinuous, containing a strong sh(x^k, or group of strains propagating at the same 
velocity. 

Another type of complication arises when the string is allowed to undergo arbitrary spacial 
motion. Even though the theory is restricted to one-dimensional strain, the velocity at every point 
of the string is a three-dimensional vector which can be resolved into tangential and normal com- 
ponents representing longitudinal and transverse motions respectively. Depending on the ini- 
tial shape of the string and the boundary conditions, the two wave motions may be coupled in 



' N.K.C.-N.B.S. Postdccloral K<-si(lftit Hcscarcli Associate. 1966-68. 
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the sense that no simple wave solution for each motion can be found. Two characteristic veloci- 
ties are involved: different values of strain in the longitudinal wave are propagated at velocity 
c given by (1), and changes in configuration of the string are propagated at different values of the 
velocity P given by 



For a detailed discussion of this rate-independent theory as applied to string motion, see Cristescu 
[1, 2V 

The comphcation in the coupling of the two motions is avoided when one considers an infinitely 
long string, initially straight and subjected to a constant-velocity transverse impact. Simple wave 
solutions were reported by Smith, McCrackin, and Schiefer [3], and were recently given in another 
form by Schultz, Tuschak, and Vicario [4]. Attempts to interpret the results of transverse impact 
experiments on materials such as nylon and polyester yarns, and strips of natural rubber, based 
on the rate-independent theory, however, were not completely satisfactory. In particular, two 
basic predictions of the theory were not verified: 

(a) The quasistatic and dynamic stress-strain curves of the nylon and polyester samples 
tested by Petterson and Stewart [5, 6J and by Smith, Fenstermaker, and Shouse [7J, have portions 
concave towards the tension axis, implying that strong shocks should exist in the strain distribu- 
tions. Neither group noted the presence of any strong discontinuities in their experimental data. 

(b) For each impact velocity, the theory yields a unique relation between the strain and the 
longitudinal wave velocity c(e). It is an inherent feature of the rate-independent theory that this 
unique relation is valid for all impact velocities. Experimental data failed to support this assertion 
after appropriate consideration was given to measurement uncertainties. For additional discussion, 
see Smith and Fenstermaker [8]. 

In order to gain insight for interpreting these effects it is useful to study the behavior under 
impact conditions of a simple linear viscoelastic model, such as one consisting of a spring in 
parallel with a Maxwell element.*^ This model has already been used by Morrison [9] and Smith [10] 
to study wave propagation in filaments impacted longitudinally. By analogy with results obtained 
in these theoretical studies, the effects of viscoelasticity on experimentally obtained strain, tension, 
and particle velocity distributions are better understood. By extending the use of this model to 
the transverse — impact problem, further enlightenment will result. In a viscoelastic string it is 
not possible to uncouple the mutual influence of the two types of wave as was done in the rate- 
independent theory under some restrictive assumptions. Thus one must study theoretically such 
additional effects of viscoelasticity as changes in the transverse wave profile, changes in the trans- 
verse wave propagation velocity, and secondary changes in the strain, tension, and particle velocity 
distributions. 

A detailed analysis of such a problem with a proposed numerical scheme suitable for the 
calculation of strong discontinuities was pubHshed by Smith [11]. Since pubhcation several errors 
have been found with the consequence that the proposed numerical scheme is not vaHd, although 
the basic formulation of the problem and the mathematical method of solution remain apphcable. 
In this paper these errors are corrected and the numerical scheme revised. In addition the problem 
is discussed from a different viewpoint and an alternate numerical scheme presented. Illustrative 
numerical calculations are also given. 



- Fijiures in brackets indicate the literature references at the end of this paper. 

•' It should be noted that for three-dimensional problems there are objections to the use of constitutive equations characterizing the behavior of spring and dashpot 
models. Unless properly formulated, such equations, containing time derivatives of the stress and strain, violate the so-called "principle of material indifference" 
as stated by Truesdell and Toupin |12|. However, the restrictions imposed by this principle are automatically satisfied in the one-dimensional problem considered 
here. 
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2. Notation 

The symbols used are the same as those used in reference [11] except for the symbol K repre- 
senting the spring constant. 

m initial density of the string in mass per unit length. 

K spring constant of the string in force units. 

X dimensionless parameter having a value hetvs^een and 1 such that XK is the spring constant 

of the Maxwell element constituting one branch of the viscoelastic model. 
T relaxation time of the Maxwell element. 

c velocity of the longitudinal wavefront as given by eq (1); for the three-element model it is 

given hy c=v Kim. 
T tension in the string in force units. 

P velocity related to the propagation of transverse waves as given by eq (2). 

V transverse impact velocity. 
t time coordinate. 

X Lagrangian space coordinate fixed to the string. 

x\ y laboratory or observer's coordinates to be referred to as horizontal and vertical coordinates, 

respectively; the x and x coordinates of the string coincide in the initial unstrained state. 
^, 17 horizontal and vertical components of displacement of a point x on the string relative to its 

original position: 

(3) 

X =x^-^[x. t): y =r)(x, t). 

u, V horizontal and vertical coVnponents of velocity in x , y' space of a point x on the string: 

— ^— M. _(}y__'^ 

dt dt' in ()t 

a, jQ first partial derivatives of the laboratory coordinates of a point x on the string with respect 
to the Lagrangian coordinate x: 

^ = ^=1 + ^,^ = 1^ = ^. (5) 

c)X OX dx c)X 

ds infinitesimal length element of the deformed string: 



ds = \/{dx')-^{dy')-=Va-^p-dx. (6) 

€ Lagrangian strain at x on the string 

ds — dx . / ., , ^., , 
^ = —] =Va- + i8--l. (7) 

6 angle that each infinitesimal length element makes with the horizontal x' direction: 

„ dx' dx a . .. B 
cos = -—-p = j-—; sme = j^' (8 

c)x ds 1 + e 1 -h e 

w flow or particle velocity of a point x on the string with respect to an inextensible coordinate 

coinciding with the string and fixed to it at points in advance of the longitudinal-wave front; 
that is, the rate of extension of the portion of string between point x and the longitudinal- 
wave front. This velocity satisfies the differential equation: 
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dw = du cos 6-\-dv sin 6 = 



1 + e 
3. Formulation 



du- 



1 + e 



dv. 



(9) 



The viscoelastic model consists of a spring of spring constant (1— \)/C in parallel with a 
Maxwell element of spring constant \K and relaxation time r. If the quantity mc- is substituted 
for K, the constitutive equation for this model takes the form 



., de . {l — k)mc- c)T . I ^ 

dt T dt T 



For plane motion the following two equilibrium equations apply: 

aT 



(10) 






dt c)x ax 



1-he 
f3T 



1-he 



The equations of continuity are 



du 

dx~ 


da 
~ dt 


dx 


M 

dt 



(11) 



(12) 



(13) 



(14) 



We shall consider only the waves propagating in the positives direction. The solution therefore 
of the set of eqs (7), (9), and (10) through (14) must satisfy the initial conditions: for f=0, x > 



u = v = w = li = €= T=0,a= 1, 

and the boundary conditions: for t> 0, x = 

u = 0, v=V, 

and ior t> 0, x-*-^ 

u = v = w — P = e=T=0, a= 1. 



(15) 



(16) 



(17) 



In addition to the initial and boundary conditions, certain compatibility conditions also apply. 
At the longitudinal wave front e, w, and T are discontinuous and must satisfy the jump relations 



[w\ = -c[e] 
[T] = — mc[w] 



(18) 
(19) 



where the bracketed quantities represent jumps in value occasioned by the passage of the wave 
front at the point under consideration. At the transverse wave front u, v, a, /3, are discontinuous 
and must satisfy the jump relations 

[u] = -P[a-l] (20) 



[v] = -P[(il 
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(21) 



4. Explicit Method of Solution 

The variables /^ v, a, (3 can be eliminated from eqs (7), (9), and (11) tbrouj^li (14) to obtain the 
equations 

aw 1 aT 

(22) 
(23) 



dw 
dt ~ 


m 


dT 

dx 


dw\ 
dx " 







These equations and the constitutive eq (10) form a set having characteristics given by 

dx 



(It 



(24) 



and 

dx 

dt 



= ^' i:2S) 



Thus, the strain e, tension T, and particle velocity w can have discontinuities only along the charac- 
teristic curves jc=^±cf-f a, and x=b, in which a and h are constants of integration. In particular, 
e, r, and w are discontinuous at the longitudinal wave front, but are continuous at the transverse 
wave front. 

The following efjuations hold along the characteristics: For the ''+c'' characteristic 

mcdw — dT = [T — {I — X)mc^e] — ; (26) 



for the "— c" characteristic 



for the "x = b'^ characteristic 



fncdw-\- dT = — [T — (l — X)mc-e] — ; (27) 

T 



mc'de - dT= [7- (1 - X)mc'€] — • (28) 

T 



Equation (28) is the same as the constitutive eq (10). 

The values of the strain, tension, and particle velocity at the longitudinal wave front can be 
found by solving the set of equations (18), (19), and (26). The results are 

e(cr, = ^oe-^^^/2r) ^29) 

T(ct, t) = mc'€oe-^^fl^''> (30) 

w(ct, t) =-ceoe-<^^/2^\ (31) 

where eo is the strain 6(0, 0) resulting at the instant of impact. To find the vahie of eo, the com- 
patibility conditions (20) and (21) are expressed in tlu^ form 

u = w^P{\^€-a) (32) 
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v = -Pp (33) 

in which the values of u, v, a, and j8 are understood to be the values obtained just after passage of 
the transverse wave front, and the continuous quantities P, e, and w are also evaluated at the 
wave front. Applying the boundary conditions (16) to eqs (18), (32), and (33), and using eq (7), one 
finds that at the point of impact the relation V^ = Pl(\ -\- e^Y — [(1 + e())Po — ceo]" must be satisfied. 
From eqs (2), (18), and (19), one obtains Po = cVeo/(l + eo) which, when substituted into the rela- 
tion above, gives 



P/c^ = 6o(l + €o)-[V€o(l + 6o)-eo]^ (34) 

The value of eo is obtained by solving eq (34). 

If velocity P is defined as in eq (2), it is obvious that eqs (11) through (14) form a set involving P 
and the dependent variables a, r, a, and f^. This set of equations has characteristics given by 

f=^^ '^^^ 

and along these characteristics the following relations hold: For the '' + P" characteristic, 

f3{du-Pda)-a{dv-Pdp) = 0, (36) 

and for the "— P" characteristic, 

fBidu-\-Pda)-a{dv-\-Pdf3)=0. (37) 

The solution of the system of equations (7), (9), (26), (27), (28), (36), and (37), is found by numeri- 
cal integration along the set of characteristics plotted in the Lagrange diagram, figure 1. In this 
diagram the straight line characteristics are solutions of the equation (ix/G?^ = ±c; the characteristic 
passing through the origin, and having the equation x = c^ corresponds to the discontinuity marking 
the wave front of the strain, tension, and particle velocity distributions. Thus, in region between 
this characteristic and the x axis the strain e, tension T, and particle velocity w are zero. Along the 
region 1 side of the characteristic the values of e, T, and w are given by eqs (29), (30), and (31), 
which may be thought of as boundary conditions. 

The curved line passing through the origin represents a solution of the equation dxldt = P. 
The exact course of this characteristic is not known in advance, as P is a function of e and T. This 
characteristic corresponds to the front of the wave of transverse motion that propagates along the 
filament; thus, in region 1 between this characteristic and the characteristic ::t; = ci the velocity u 
is equivalent to the velocity w, the variable a is equal to 1 + e, and the velocity f and variable 13 
are equal to zero. The jump relations (32) and (33) holding along the transverse wave front constitute 
a second set of boundary conditions. A third set of boundary conditions expressed by eq (16) holds 
along the time axis. 

Calculation of values for u, f , w, a, /3, €, T at the mesh points i;: the Lagrangian diagram (using 
fig. 1 as an example) proceeds systematically in the following steps: 

1. Calculate €o using eq (34). 

2. Calculate w, e, T for mesh points along the characteristic x = ct, using eqs (29), (30), and 
(31). 

3. Calculate w, a, j8, €, T at point 1, and u, v, w, a, /3, e, T at point 2 using a subroutine to be 
described later. 

4. Calculate u, v, w, a, /3, e, T at successive points 3, 4, 5, 6. 

5. Calculate w, e, T at successive points a, 6, c, etc. 

6. Continue the calculation as indicated. 
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F IGURE 2. Mes/i situations encountered during integra- 
FlGURE 1. Characteristics network used in the explicit tion along the characteristics plotted in figure 1. 

method of solution. 

The straight Hne characteristics of positive and nejiative slope are solutions 
of th«- eq dx/dt = ±c. The characteristic separating region 1 fnnn region 
gives the location of the longitudinal wave front. The curve characteristic 
separating region 2 from region 1 and satisfying the eq dxldt = l' gives the 
location of the transverse wave front. 

Five mesh situations depicted in figure 2 are encountered in the numerical integration process. 
The simplest of these is that of figure 2c. Here, points 1 and 5 lie at the intersections of c charac- 
teristics, and the characteristic between points 2 and 4 may be either a c or a P characteristic; 
line (3, 5) is a constructed vertical line. Values of it;, e, and T at point 5 are found from data known 
at points 1, 2, and 4 by solving eqs (26), (27), and (28) written respectively in the form 



mew— T— mcw\ — T\^ [T\:^ — (1 — \)mc-E\-^\{t:, — tij/r 
mew +7*— mew\ -\-Ta — [T^:, — (1 — \)me-E47,\{h — tA)lT 
me'-e — T— me'e-A — T-.i -h fr^r, — ( 1 — \)mc-E:ir^\{t:, — t:i)lT. 



(38) 
(39) 
(40) 



in these equations the symbol T\:,, for example, represents the average value of {T\-\-T)l2. In the 
first trial calculation T\:, is set equal to T\\ the calculation is repeated twice using improved values 

ofr„. 

Figure 2a depicts a situation found in meshes low on the t axis. In the lowest mesh the trans- 
verse wave front characteristic passes through the origin so that points 1 and 2 coincide, but at 
later times the more general form depicted may occur. Values of u, v, w, a, jS, e, T, jc, and t at point 

6 and w^ a, j8, e, and T at point 7 are required. The large set of equations involved is most con- 
veniently solved by an iterative process in which it is assumed initially that the value of z^ at point 

7 is equal to the value of w at point 1. This enables e, and T at point 7 to be calculated using eq 
(28) along the time axis (1, 7) and eq (27) along characteristic (4, 7) (segment (6, 7) in subsequent 
iterations). 

The transverse wave front characteristic is constructed with slope 1/P evaluated from point 2 
data. The vertical line (3, 6) and the characteristic (5, 6) are then constructed and values of z^, e, and 
T at points 3 and 5 determined by interpolation. Values of w, e, and T are calculated at point 6 
using eq (26) along (5, 6), eq (27) along (4, 6), and eq (28) along (3, 6). Values of w, i;, a, and (i are found 
using eq (36) along the transverse wave front characteristic (2, 6) and the compatibility relations 
(32) and (33) plus eq (7) applied at point 6. Characteristic (7, 8) is then constructed with slope 
(— 1/P) evaluated from point 7 data. The values of a and ^ at point 7 are found l)y applying (U[ (37) 
along (7, 8) and eq (7) at point 7. The iteration cycle is completed by findings; at point 7 through use 
of eq (9) along line segment (6, 7). Improved values of the required quantities are obtained by 
iterations which are continued until successive values of w at point 7 are found to differ by a suffi- 
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ciently small amount. During the calculation subiterations of the type discussed for situation 2c are 
also carried out. 

Figure 2b depicts the situation encountered in tracing out the course of the transverse wave 
front characteristic; values of //, v, w\ a, /3, e, T, x. and t are required at i)oint 7. In this case values 
of w, e, and T are found from relations holding along (6, 7) and the two constructed lines (2, 7) and 
(5, 7). Values of a, v, a, and /S are found using eq (36) along (3, 7) and applying the compatibility 
relations (32) and (33) plus eq (7) at point 7. 

Situation 2d occurs when values oi u, v, a, and 13 at point 6 are calculated from data known at 
points 1, 3, and 5, and from values of w, e, and T at point 6. The characteristic between points 
3 and 5 may be either a c characteristic or the transverse wave front characteristic. In this case, 
eqs (36) and (37) are applied respectively along the constructed P characteristics (2, 6) and (4, 6), 
eq (9) is applied along the c characteristic (5, 6), and eq (7) is apphed at point 6. 

In situation 2e values of w, a, /3, €, and T at point 4 on the time axis are calculated from data 
at points 1 and 3 on a c characteristic. In the solution of this problem, eq (28) is applied along the 
time axis, eq (37) along the constructed P characteristic (2, 4), eqs (9) and (27) along the c charac- 
teristic (3, 4), and eq (7) at point 4. 

In order to find the configuration of the filament in laboratory {x\ y') coordinates, it is desir- 
able to know the values of ^ and r] at each of the mesh points in the Lagrange diagram. Along the 
longitudinal wave front the values of ^ and r] are zero, and along the transverse wave front the value 
of 7} is zero. The value of ^ is continuous across the transverse wave front. In region 1 the value 
of 7} is zero and the value of ^ is found in terms ofw and e by integrating along a negative c charac- 
teristic using the relation 

dt cU dt dx 
which becomes, after appropriate substitutions 

d^={w-ce)dt, (42) 

In region 2, ^ and tj are found by integrating along a negative c characteristic using the relations 

d^=[u-c{a-l)]dt (43) 

d'r] = {u-c^)dt. (44) 



5. Errors in a Previous Soiution 

The scheme just presented has been outlined in more than the usual detail because it avoids 
several errors occurring in a previous solution [11]. Most of these errors were due to improper use 
of eq (9). It should be noted that eq (9), upon substitution of the boundary values for u and v at 
x=0, reduces to the relation {dwldt)x=o= 0. This implies erroneously that w is constant along the 
time axis. The mistake occurs because in this appHcation of eq (9) w has been expressed relative 
to a coordinate system fixed to the string at the point of impact, but the value of w desired is that 
relative to a system fixed to the string at points in advance of the strain wave front. A similar error 
results if eq (9) is applied along the transverse wave front characteristic. 

Another erroneous result is obtained if, in situation 2e, eq (9) is applied along the negative P 
characteristic to find the value of w at point 4. In this example the value of w at point 2 is obtained 
by interpolation between values at points 1 and 3. Thus, the relation (9) is split into two components 
one of which is appHed along the t axis, situation 2a provides a similar opportunity to make this 
mistake. 
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Equation (7) can be expressed in the differential form 

ada + /3c/^ = ( IH- e We (45) 

which is more convenient for calculation. Its use in this form, however, is unnecessary, and pro- 
duces a small cumulative error. 

If eqs (36) and (37) are combined with eqs (9) and (45), the following relations are obtained 
along the "'-{' P'' characteristics: 

du - Pda = -^ ( dw — Pde) (46) 

dv - Pdp = -^ {dw - Pde); (47) 

and along the "— P" characteristics: 

du -f Pda = ——- (dw + Pde ) (48) 

dv -h Pdp = y|- idw + Pde). (49) 

However, these relations, used previously [11], have yielded inaccurate and often erroneous results. 
They are, for instance, invalid along the transverse wave front characteristic when used in connec- 
tion with eq (7) or (45), and invalid along the ' — P" characteristic in situations 2ci and 2e. In situ- 
ation 2d eqs (46) through (49) have yielded inaccurate results, especially in those meshes involving 
the transverse wave front characteristic. 

6. Coupling Between Longitudinal and Transverse Waves 

It is interesting to note that of the two wave phenomena present, the longitudinal wave is 
described in terms of the quantities 6, ii\ and T, and the transverse wave in terms of the ([uantities 
//, r, a, and /3. It has been shown that the differential equations describing the problem can be 
reduced to a set involving ii\ e, and 7' only; thus, it would seem that values of w, 6, and 7' could 
be found separately, and in a second step values of u, v, a, and ^ found in terms of them. Un- 
fortunately, however, the boundary conditions are expressed in terms of the variables a and v 
so that in general this decouphng is impossible."* 

From a physical point of view the effects of coupling between the variables describing the 
the longitudinal and transverse waves can be described as follows: Changes in the strain, tension, 
and particle velocity distributions due to viscoelasticity induce changes in the configuration and 
propagation velocity of the transverse wave which react back to modify the strain, tension, and 
particle velocity distributions. Thus, the viscoelastic effects can be considered as both primary 
and secondary. 

The distributions resulting from the primary effect, are the same as those resulting from a 
longitudinal impact at constant velocity, in which the strain increases, the tension decreases, and 
the particle velocity stays constant at the point of impact. In the secondary or coupling effect the 
particle velocity at the point of impact changes with time, and associated changes occur in the 
strain, tension, and particle velocity distributions. If viscoelasticity is absent, as when the constitu- 
tive equation is of form T = T{e), and if transverse impact occurs at constant velocity, only the 
primary effect is involved, and no changes occur in the longitudinal and transverse waves. In this 
case the problem may be considered as effectively decoupled. 



Mm iIk- previous papi-i |ll| <■(] (<>) was applied aloiij: llic / axis and tin- emmcniis result k =(i constant used as a boundary condition t(. achieve deeouplir 

209 



7. Implicit Method of Solution 

In the method of solution discussed above, expHcit systems of equations are used to evaluate 
the unknowns at successive points in a characteristics network. Alternatively, values of the un- 
knowns at points, along the string can be calculated simultaneously by the following algorithm, 
based on an implicit finite-difference formulation for nonlinear wave propagation problems [13]. 

In order to simplify the calculation, the mathematical problem is reformulated in terms of the 
five unknowns u, v, a, /3, and T. The two additional unknowns, e and w, are defined by eqs (7) 
and (9), and can be found separately. The equations are the same as those used in the first method; 
namely, eqs (26) and (27) holding along the ±c characteristics, the constitutive eq (28), and eqs 
(36) and (37) holding along the ±P characteristics, except that the variables e and w have been 
evaluated in terms of the variables u, v, a, and j8, and do not appear. 

The mesh system used is depicted in figure 3. As before, there are three regions of interest: 
region for points in advance of the longitudinal wave front, region 1 for the portion of the longi- 
tudinal wave in advance of the transverse wave front, and region 2 for points within the transverse 
wave. The differential equations have different forms in regions 1 and 2. In region 1 there are no 
P characteristics and the variables v and (3 are both equal to zero; thus only eqs (26), (27), and (28) 
with a— I substituted for e and u substituted for w are used here. The more general forms of the 
equations with Va- 4-/3^ — 1 substituted for e and ( adu -\- I3dv)l\^d^ + /3^ substituted for dw are 
used in region 2. Along the longitudinal wave front the values of u, a, and T are found from eqs (29), 
(30), and (31) by substituting a— 1 for e and u for w, and along the t axis the values oi u and v are 
given by the boundary condition (16). Along the transverse wave front characteristic four com- 
patibility conditions are needed. These are: 



U2 = U\—P{a2 — ai) 

V2 = -PIB2 
T2 = T, 



(50) 

(51) 
(52) 



«i + iSi = «? 



(53) 




X 

< 



X AXIS 

Figure 3. Mesh system used in the implicit method of 
solution. 

The Ax increments of the mesh are constant in rejiion 2. The Ax increments 
in region 1 and A/ i^crement!^ in refiions 1 and 2 are variable, dependinti 
on the value of the velocity /^ at the transverse wave fmnt. 
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where the subscripts refer to values of the variables in the regions 1 and 2. Equations (50) and (51) 
are the same as the compatibility conditions (20) and (21), and eqs (52) and (53) merely state mathe- 
matically that the tension and strain are continuous across the transverse wave front as was 
shown earlier in the text. Alternatively, eqs (52) and (53) may be obtained by equating horizontal 
and vertical components of the force across the transverse wavefront [3] or by the matrix methods 
expounded by Jeffrey and Taniuti [14]. 

In the mesh system the distance increment Ax is maintained constant in region 2, and a variable 
time increment is defined as ^t=^xlP where P is the transverse wave front velocity averaged 
between points (x, t) and (x + Ax, t-\- At). In region 1 the distance increment depends on the value 
of P and thus is variable also. The value of P as a function of x and t is not known in advance, so 
the mesh system can only be constructed in stages as the calculation proceeds. 

The calculation is carried out in a series of timewise steps in which the values of w, f , a, j8, 
and T at each mesh point corresponding to t-\- At, are found in terms of values of these quantities 
known at each mesh point corresponding to time t. In the calculation each of the timewise steps is 
iterated several times in order that the trajectory of the transverse wave front (and value of P) be 
known with sufficient accuracy. 

Two types of calculation, represented schematically in figure 4, occur in each step. The sim- 
plest of these, figure 4a, is used in region 1 to calculate values of //, a, and T at point 6 in terms of 
values known at points 1, 2, and 5. In this calculation vertical line (3, 6) and fine (4, 6) of slope 
( — 1/c) are constructed, and the values of //, a, and T at points 3 and 4 determined by interpolation. 
The required values at point 6 are found by solving the equations holding along the characteristics 
(1, 6), (3, 6), and (4, 6). The calculation is applied sequentially starting at the longitudinal wave front 
characteristic and working towards the transverse wave front characteristic. 

The subcalculation, figure 4b, is used to find values of Wi, U), V\, v-i, ol\, ol^, and T at point a 
on the transverse wave front and values of //, a, and T at the mesh points in region 2 corresponding 
to time t-\-At. The subscripts refer to regions 1 and 2 separated by the transverse wave front 
characteristic. For the seven unknown quantities at point a there are seven corresponding equa- 
tions, four holding along the characteristics (/>, a), (c, a), (J, a), and (e, a) plus the three compati- 
bility conditions (50), (51), and (53). Values of a, v, a, /3, and T are not known at point e however, 
so the values required at point a cannot be determined explicitly by solving these seven equations. 
To find the five values required at point e it is necessary to know the corresponding values at point/, 
and this requires the introduction of five equations holding along the c characteristics (^, /) and 
(7,/), the P characteristics (h,f) and (i,/), and the fine {d,J). Similarly, additional equations are 
required for each of the mesh points in region 2. The final point p lies on the t axis and is evaluated 
from three equations holding along (^, p), (r, /;), and (s, p) plus the boundary conditions m = 0, and 
f = V. In this way a large system of simultaneous equations is obtained, the solution of which gives 
the required values at each of the mesh points. 








Fl(;URE 4. Mesh situations encountered in the implicit 
method of solution 
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This alternate formulation has several advantages. A fundamental set of equations is used in 
which the number of unknowns is equal to the number of relations holding along the characteristics. 
This set of five equations in u, v^ a, /3, and T is well posed for the prescribed initial and boundary 
conditions. Moreover, implicit schemes of this type are known to be stable at large values of time. 
Ordinarily, implicit schemes do not require small time steps in regions where the characteristics 
are nonlinear, but in this application the presence of a moving boundary of unknown curvature 
(transverse wave front characteristic) may limit this latter advantage. 

8. Calculated Results 

The explicit method has been used to calculate the behavior of the viscoelastic string for 
various values of the parameters X, m, c, and r characterizing the model and for various values of 
the transverse impact velocity V. The results for one of these calculations with parameters X= 0.2, 
m= 100 gm~^ T= 0.002 s, c= 100 ms~^ and V= 50 ms"' are presented here. The values chosen for 
the parameters roughly simulate the material properties of an elastomer such as rubber when sub- 
jected to transverse impact. 

The distributions of the particle velocity, strain, and tension calculated in this example are 
shown in figures 5, 6, and 7, respectively. The ordinate in figure 5 is a dimensionless quantity, the 
ratio of the particle velocity at a Lagrangian distance x along the string to the particle velocity at 
the instant of impact, w/wo. The distribution of particle velocity in the wave is shown by the solid 
line curves for times of 1, 2, 3, 4, and 5 ms after impact. The dashed lines indicate the velocity 
distribution for the case of no relaxation {k=0). The distributions are similar to those resulting 
from longitudinal impact at constant velocity [10] except that in this example the particle velocity 
at the origin increases with time sHghtly, whereas in the case of longitudinal impact it necessarily 
remains constant. This illustrates the effect of coupling between the longitudinal and transverse 
waves when viscoelasticity is present. 

The dimensionless ordinate in figure 6 is the ratio of the strain at a distance x to the strain at 
the instant of impact, e/eo. The distributions depicted are similar to those resulting from longitudinal 
impact at constant velocity. The strain at the point of impact increases with time and approaches a 
limiting value, but this limiting value is greater than the value of eo(l — X)"^^"= l.HSeo which is 
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Figure 5. Distribution of particle or flow-velocity along the string for times of 
1,2,3, 4, and 5 X 10~^ s after transverse impact. 

Calculation parameters are X = 0.2, m= 100 gm~', c= 100 ms"', t = 0.002 s, and ^' = 50 ms~'. Dash lines give 
distributions for an elastic string (\ = 0). The particle velocity is expressed as the ratio of the velocity it at distance 
X to the velocity uo attained at the instant of impact. 
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Figure 6. Distribution of strain along the string for times of 1, 2, 3, 4, and 5 X 10 "^ s 
after transverse impact. 

Calculation parameters are \ = ().2, m^ 1(K) fim-'. c= 100 ms"', t = 0.002 .s, and F = 50 ms-'. Dash lines give 
distributions for an elastic string (X = 0). The strain is expressed as the ratio of the strain e at distance x to the 
strain co attained at the instant <tf impact. 
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Figure 7. Distribution of tension along the string for times of 1,2, 3, 4, and 5 X 10^'^ s 
after transverse impact. 

Calculation parameters are \ = ().2, m = 100 gm^', c=100 ms"', t = 0.002 s, and V = 50 ms '. Dash lines give 
distributions for an elastic string (A = 0). The tension is expressed as the ratio of the tension Tat distance x to the 
tension T» attained at the instant of impact. 



attained in the case of longitudinal impact [10]. The tension distributions at various times after impact 
are depicted in figure 7. These distributions also are similar to those resulting from longitudinal 
impact at constant velocity. The tension at the point of impact decreases with time and approaches 
a hmiting value which is less than the value of 7^(1 — \)'/- = 0.894r() attained in the case of longi- 
tudinal impact. 
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The values of particle velocity, strain, and tension attained at the instant of impact are — 27.50 
ms~S 0.2750, and 2.750 X 10- N. The initial value of the Lagrangian transverse wave front velocity 
is 46.44 ms~^ but after 5 X 10"'^s it slows to 44.00 ms ^ To an observer in the laboratory U', y' sys- 
tem) the initial velocity is (1 +6())P() + i(;() = 31.71 ms~^ and after 5X lO'^s the wave front has trav- 
eled a distance x' = x-\- ^= 0.1525 m. If the wave front had traveled at the initial velocity, it would 
have gone 0.1586 m. At the instant of impact the slope of the transverse wave, tan 6 = (B/a, is equal 
to — 1.577, but after 5X 10 -^s it decreases to —1.646 at the point of impact and —1.627 at the 
transverse wave front. Thus, the wave slows down as it propagates and gets shghtly steeper and 
becomes slightly curved. However, for this calculated example these secondary coupling effects 
are small and would not be easily discerned by an observer in the laboratory. 

9. Summary and Conclusions 

This paper has described two methods for calculating wave propagation effects in a linearly 
viscoelastic string subjected to transverse impact at constant velocity. Results have been presented 
for a calculated example. 

Viscoelasticity exerts both primary and secondary effects on the distributions of stress, strain, 
and particle or flow velocity along the string that result from constant-velocity transverse impact. 
The distributions resulting from the primary effect are the same as those resulting from a longitu- 
dinal impact at constant velocity, in which the strain increases, the tension decreases, and the 
particle velocity stays constant at the point of impact. In the secondary or coupling effect visco- 
elastic relaxation causes the particle velocity at the point of impact to increase with time, and 
associated changes to occur in the strain, tension, and particle velocity distributions. 

The secondary effect also causes the transverse wave front to propagate at a gradually decreas- 
ing velocity. In addition, the transverse wave profile gets steeper and becomes gradually curved. 
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